program test2x2
	
	use QMCTEST
	
	implicit NONE
	
	integer(4), parameter		:: nmax = 100000
	
	real(8)						:: beta
	integer(4)					:: H(L)
	integer(4)					:: lmc, nmc
	integer(4)					:: ip
	real(8)						:: sz
	integer(4)					:: ncnt(0:1)

	integer(4), parameter		:: ND = 16
	integer(4)					:: lmax
	integer(4)					:: nc
	real(8)						:: zsum, nsum, n2sum
	integer(4)					:: hdg(ND)
	integer(4)					:: p, q
	real(8)						:: d, zz, ee
	
	real(8)						:: evals(ND), E(ND), Z(ND), E2(ND)
	real(8)						:: Zs, Es, Eavr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	hdg = [	4, 0, 0, 0, 0,-4, 0, 0, 0, 0,-4, 0, 0, 0, 0, 4]
	
	call Init()
	
	evals = jj * hdg + cj
	
	lmax = 2**L
	zsum = 0.0d0
	
	beta = 1.0d0
	
	do p = 1, 16
		
		zz = 0.0d0
		ee = 0.0d0
			
		do nmc = 1, lmax
			
			lmc = nmc - 1
			do ip = 1, L
				H(ip) = mod(lmc,2)
				lmc = lmc / 2
			end do
			nc = count( H(1:L) == 1)
			d = 1.0d0
			do q = 1, L
				if ( H(q) == 1) d = d * evals(p)
			end do
			
			zz = zz + d * fact(nc) * (-beta)**nc
			ee = ee + d * fact(nc) * (-beta)**nc * nc
			
		end do
		
		Z(p) = zz
		E(p) = ee
		
	end do
	
	continue
	
	Zs = sum( Z(:))
	Es = sum( E(:))
	
	Eavr = - Es / Zs / beta - cj
!	Eavr = Es / Zs - cj
	
	stop
	
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	beta = 1.0d0
	H = 0
	ncnt = 0
	call Init()
	do nmc = 1, nmax
		call RANDOM_NUMBER( sz)
		if ( sz > 0.5d0) then
			call LocalDuoStep(beta,H)
		else
			call LocalUniStep(beta,H)
		end if
		forall( ip = 0:1) ncnt(ip) = ncnt(ip) + count( H(1:L) == ip)
	end do
	
	continue
	
end program test2x2